1 Introduction

Reviewer #2 requested the following analysis for the T-cell activation and culturing experiment:

  • Characterization of the transcriptomic changes associated with exposure to anti-CD3 antibody.
  • Compositional analysis of cell types with and without anti-CD3 antibody treatment.

2 Pre-processing

2.1 Package loading

library(scater)
library(Seurat)
library(ggpubr)
library(purrr)
library(kBET)
library(viridis)
library(ggmap)
library(cowplot)
library(tidyverse)

2.2 Source script with functions

source("bin/utils.R")

2.3 Load data

t_act_male1 <- readRDS("../3-T_cell_activation/results/R_objects/t_act_Seurat_male1.rds")
t_act_rep2_l <- readRDS("../3-T_cell_activation/results/R_objects/t_act_Seurat_list_reannotated_rep2.rds")

2.4 Merge seurat objects and homogenize annotations

# Donor 1
t_act_male1$`0M`$cell_type <- factor(
  t_act_male1$`0M`$cell_type,
  levels = c("CD4 T", "Cytotoxic", "Monocyte", "B", "FCGR3A Monocyte", "Dendritic Cell")
)
levels(t_act_male1$`0M`$cell_type) <- c("CD4 T-cell", "Cytotoxic", "CD14 Monocyte", "B-cell",
                              "FCGR3A Monocyte", "Dendritic Cell")
t_act_male1$`2M`$cell_type <- factor(
  t_act_male1$`2M`$cell_type,
  levels = c("CD4 T", "Cycling", "Cytotoxic", "B")
)
levels(t_act_male1$`2M`$cell_type) <- c("Activated CD4 T-cell", "Cycling CD4 T-cell",
                                        "Cytotoxic", "B-cell")
t_act_male1$`0M`$is_cultured <- rep("uncultured", ncol(t_act_male1$`0M`))
t_act_male1$`2M`$is_cultured <- rep("cultured", ncol(t_act_male1$`2M`))
t_act_male1 <- merge(
  x = t_act_male1$`0M`,
  y = t_act_male1$`2M`,
  add.cell.ids = c("uncultured", "cultured")
)
t_act_male1$donor <- rep("donor1", ncol(t_act_male1))

# Donor 2
t_act_rep2_l$`0_rep2_F2`$cell_type <- factor(
  t_act_rep2_l$`0_rep2_F2`$cell_type,
  levels = c("CD4 T-cell", "CD8 T-cell", "NK", "B-cell", "Monocyte", "FCGR3A Monocyte", "Unknown")
)
levels(t_act_rep2_l$`0_rep2_F2`$cell_type) <- c("CD4 T-cell", "CD8 T-cell", "NK",
                                                "B-cell", "CD14 Monocyte", "FCGR3A Monocyte",
                                                "Unknown")
t_act_rep2_l$`1_rep2_F2`$cell_type <- factor(
  t_act_rep2_l$`1_rep2_F2`$cell_type,
  levels = c("Activated CD4 T-cell", "Cycling CD4 T-cell", "CD8 T-cell", "NK", "B-cell", "Unknown")
)

t_act_rep2_l$`0_rep2_F2`$is_cultured <- rep("uncultured", ncol(t_act_rep2_l$`0_rep2_F2`))
t_act_rep2_l$`1_rep2_F2`$is_cultured <- rep("cultured", ncol(t_act_rep2_l$`1_rep2_F2`))
t_act_female2 <- merge(
  x = t_act_rep2_l$`0_rep2_F2`,
  y = t_act_rep2_l$`1_rep2_F2`,
  add.cell.ids = c("uncultured", "cultured")
)
t_act_female2$donor <- rep("donor2", ncol(t_act_female2))

# Donor 3
t_act_rep2_l$`0_rep2_F3`$cell_type <- factor(
  t_act_rep2_l$`0_rep2_F3`$cell_type,
  levels = c("CD4 T-cell", "CD8 T-cell", "NK", "B-cell", "Monocyte", "FCGR3A Monocyte")
)
levels(t_act_rep2_l$`0_rep2_F3`$cell_type) <- c("CD4 T-cell", "CD8 T-cell", "NK",
                                                "B-cell", "CD14 Monocyte", "FCGR3A Monocyte")
t_act_rep2_l$`1_rep2_F3`$cell_type <- factor(
  t_act_rep2_l$`1_rep2_F3`$cell_type,
  levels =  c("Activated CD4 T-cell", "Cycling CD4 T-cell", "CD8 T-cell", "NK", "B-cell", "Unknown")
)
t_act_rep2_l$`0_rep2_F3`$is_cultured <- rep("uncultured", ncol(t_act_rep2_l$`0_rep2_F3`))
t_act_rep2_l$`1_rep2_F3`$is_cultured <- rep("cultured", ncol(t_act_rep2_l$`1_rep2_F3`))
t_act_female3 <- merge(
  x = t_act_rep2_l$`0_rep2_F3`,
  y = t_act_rep2_l$`1_rep2_F3`,
  add.cell.ids = c("uncultured", "cultured")
)
t_act_female3$donor <- rep("donor3", ncol(t_act_female3))

# Create list
t_act_l <- list(t_act_male1, t_act_female2, t_act_female3)
names(t_act_l) <- c("donor1", "donor2", "donor3")

t_act_l$donor1$cell_type <- factor(
  t_act_l$donor1$cell_type,
  levels = c("CD4 T-cell", "Activated CD4 T-cell", "Cycling CD4 T-cell", "Cytotoxic",
             "B-cell", "CD14 Monocyte", "FCGR3A Monocyte", "Dendritic Cell")
)
levels_type_female <- c("CD4 T-cell", "Activated CD4 T-cell", "Cycling CD4 T-cell",
                        "CD8 T-cell", "NK", "B-cell", "CD14 Monocyte", "FCGR3A Monocyte",
                        "Unknown")
t_act_l$donor2$cell_type <- factor(
  t_act_l$donor2$cell_type,
  levels = levels_type_female
)
t_act_l$donor3$cell_type <- factor(
  t_act_l$donor3$cell_type,
  levels = levels_type_female
)
t_act_sub_l <- purrr::map(t_act_l, function(seurat) {
  Idents(seurat) <- "time"
  seurat_sub <- subset(seurat, idents = "0h")
  seurat_sub <- seurat_sub %>%
  FindVariableFeatures() %>%
  ScaleData() %>% 
  RunPCA() %>% 
  RunTSNE(reduction = "pca", dims = 1:30) %>%
  RunUMAP(reduction = "pca", dims = 1:30)
  Idents(seurat_sub) <- "cell_type"
  seurat_sub
})

3 tSNE uncultured/cultured

donors <- c("Donor 1", "Donor 2", "Donor 3")
tsnes_is_cultured <- purrr::map2(t_act_sub_l, donors, function(seurat, titl) {
  Idents(seurat) <- "is_cultured"
  tsne <- DimPlot(
    seurat,
    reduction = "tsne",
    pt.size = 0.5,
    cols = c("#3374A1", "#E1812C")
  )
  tsne +
    labs(title = titl, x = "tSNE1", y = "tSNE2") +
    scale_color_manual("", values = c("#3374A1", "#E1812C"),
                       labels = c("Original", "Cultured")) +
    theme(plot.title = element_text(size = 13, hjust = 0.5),
          axis.ticks = element_blank(),
          axis.text = element_blank(),
          axis.title = element_blank(),
          axis.line = element_blank())
})
tsnes_is_cultured_arr <- ggarrange(
  plotlist = tsnes_is_cultured,
  ncol = 3,
  common.legend = TRUE
)
tsnes_is_cultured_arr

# Get legend
p <- tsnes_is_cultured$donor1 + theme(legend.position = "bottom")
leg <- as_ggplot(get_legend(p))
date <- Sys.Date()
# ggsave(
#   filename = str_c("../doc/figures/legends/", date, "_", "is_cultured", ".pdf"), 
#   plot = leg, 
#   width = 16, 
#   height = 5,
#   units = "cm"
# )

4 tSNE cell types

all_cell_types <- c("CD4 T-cell", "Activated CD4 T-cell", "Cycling CD4 T-cell", "Cytotoxic",
                    "CD8 T-cell", "NK", "B-cell", "CD14 Monocyte", "FCGR3A Monocyte",
                    "Dendritic Cell", "Unknown")
all_cell_types <- factor(all_cell_types, levels = all_cell_types)
palette <- c("#81324c", "#b82e57", "#e04d74", "#3a2edb", "#752bbf",
             "#c03fe4", "#bbaa2a", "#71bdd0", "green4", "hotpink2",
             "gray50")
names(palette) <- all_cell_types
tsnes_cell_type <- purrr::map(t_act_sub_l, function(seurat) {
  tsne <- DimPlot(
    seurat,
    reduction = "tsne",
    pt.size = 0.5,
    cols = palette[levels(seurat$cell_type)]
  )
  tsne +
    labs(x = "tSNE1", y = "tSNE2") +
    theme(axis.ticks = element_blank(),
          axis.text = element_blank(),
          axis.title = element_blank(),
          axis.line = element_blank())
})
tsnes_cell_type_arr <- ggarrange(
  plotlist = tsnes_cell_type,
  ncol = 3,
  common.legend = TRUE
)
tsnes_cell_type_arr

5 Cell cycle score

Garcia-Sousa I, et al. showed that the major biological process activated upon treatment with an anti-CD3 antibody is cell cycle. Let us score each cell with a cell cycle transcriptomic signature:

tsnes_s_score <- purrr::map(t_act_sub_l, function(seurat) {
  tsne <- FeaturePlot(
    seurat,
    reduction = "tsne",
    feature = "S.Score",
    pt.size = 0.5,
    cols = viridis(5)
  )
  tsne +
    labs(x = "tSNE1", y = "tSNE2") +
    theme(axis.ticks = element_blank(),
          axis.text = element_blank(),
          axis.title = element_blank(),
          axis.line = element_blank(),
          plot.title = element_blank())
})
tsnes_s_score_arr <- ggarrange(plotlist = tsnes_s_score, ncol = 3)
tsnes_s_score_arr

# Get legend
p <- tsnes_s_score$donor1
leg <- as_ggplot(get_legend(p))
# ggsave(
#   filename = str_c("../doc/figures/legends/", date, "_", "s_phase_score", ".pdf"), 
#   plot = leg, 
#   width = 16, 
#   height = 5,
#   units = "cm"
# )

6 Compositional analysis

donors <- c("Donor 1", "Donor 2", "Donor 3")
compositional_analysis_gg <- purrr::map(t_act_sub_l, function(seurat) {
  df <- seurat@meta.data %>%
    dplyr::select("is_cultured", "cell_type") %>% 
    group_by(is_cultured, cell_type) %>% 
    summarise(n_cells = n()) %>% 
    ungroup() %>% 
    group_by(is_cultured) %>% 
    mutate(pct_cells = n_cells / sum(n_cells) * 100)
  df$is_cultured <- factor(df$is_cultured, levels = c("uncultured", "cultured"))
  p <- df %>% 
    ggplot(aes(is_cultured, pct_cells, fill = cell_type)) +
      geom_col() +
      labs(x = "", y = "Percentage of cells (%)", fill = "") +
      scale_x_discrete(labels = c("Original", "Cultured")) +
      scale_fill_manual(values = palette[levels(seurat$cell_type)]) +
      theme_classic() +
      theme(plot.title = element_text(size = 13, hjust = 0.5),
            axis.title.y = element_text(size = 12),
            axis.text.x = element_text(size = 11, color = "black"))
  p
})
compositional_analysis_arr <- ggarrange(
  plotlist = compositional_analysis_gg,
  ncol = 3,
  common.legend = TRUE
)
compositional_analysis_arr

# Get legend
p1 <- compositional_analysis_gg$donor1 +
  theme(legend.position = "bottom")
p2 <- compositional_analysis_gg$donor2 +
  theme(legend.position = "bottom")
leg1 <- as_ggplot(get_legend(p1))
leg2 <- as_ggplot(get_legend(p2))
leg_list <- list(cell_types1 = leg1, cell_types2 = leg2)
date <- Sys.Date()
# walk(names(leg_list), function(leg) {
#   ggsave(
#     filename = str_c("../doc/figures/legends/", date, "_", leg, ".pdf"), 
#     plot = leg_list[[leg]], 
#     width = 16, 
#     height = 5,
#     units = "cm"
#   )
# })

7 Arrange figure

fig <- plot_grid(tsnes_is_cultured_arr, NULL, tsnes_cell_type_arr, NULL, tsnes_s_score_arr,
                 NULL, compositional_analysis_arr, nrow = 7, ncol = 1,
                 rel_heights = c(1, 0.05, 1, 0.1, 1, 0.05, 1))
# fig
# ggsave(filename = "../doc/figures/R/suppZZ.pdf", plot = fig, width = 18.5, height = 27.5, units = "cm")

8 Differential expression analysis (cultured VS uncultured)

# Merge
t_act <- merge(
  x = t_act_l$donor1,
  y = c(t_act_l$donor2, t_act_l$donor3)
)

# Process
t_act <- t_act %>%
  FindVariableFeatures() %>% 
  ScaleData() %>%
  RunPCA() %>% 
  RunTSNE(reduction = "pca", dims = 1:30) %>% 
  RunUMAP(reduction = "pca", dims = 1:30)
Idents(t_act) <- "donor"
DimPlot(t_act)

t_act$cell_type[t_act$cell_type %in% c("NK", "CD8 T-cell")] <- "Cytotoxic"
Idents(t_act) <- "cell_type"
t_act <- subset(
  t_act,
  idents = c("CD4 T-cell", "Activated CD4 T-cell", "Cytotoxic", "B-cell")
)

# DEA
t_act$cell_type[t_act$cell_type == "Activated CD4 T-cell"] <- "CD4 T-cell"
Idents(t_act) <- "cell_type"
types <- c("CD4 T-cell", "Cytotoxic", "B-cell")
t_act_types <- SplitObject(t_act, split.by = "cell_type")
dea_l <- purrr::map2(t_act_types, types, function(seurat, type) {
  Idents(seurat) <- "is_cultured"
  df <- FindMarkers(
    seurat,
    ident.1 = "cultured",
    ident.2 = "uncultured",
    test.use = "wilcox",
    logfc.threshold = 0
  )
  df
})
dea_df_l <- purrr::map2(dea_l, t_act_types, function(df, seurat) {
  mat <- as.matrix(seurat[["RNA"]]@data[rownames(df), ])
  average_expression <- rowMeans(mat)
  log2_fc <- apply(mat, 1, function(x) {
    mean_uncultured <- mean(x[seurat$is_cultured == "uncultured"]) + 0.05
    mean_cultured <- mean(x[seurat$is_cultured == "cultured"]) + 0.05
    log2(mean_cultured / mean_uncultured)
  })
  df <- df %>%
    rownames_to_column(var = "gene") %>% 
    dplyr::select("gene", "p_val", "p_val_adj") %>% 
    dplyr::mutate(average_expression = average_expression, log2_fc = log2_fc) %>% 
    dplyr::arrange(p_val_adj) %>%
    dplyr::filter(p_val_adj < 0.001)
  df <- df[, c("gene", "average_expression", "log2_fc", "p_val", "p_val_adj")]
  df
})
DT::datatable(dea_df_l$`CD4 T-cell`)
DT::datatable(dea_df_l$Cytotoxic)
DT::datatable(dea_df_l$`B-cell`)
# openxlsx::write.xlsx(dea_df_l, file = "results/tables/dea_culturedVSuncultured.xlsx")

9 Session Information

sessionInfo()
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] forcats_0.5.0               stringr_1.4.0               dplyr_0.8.4                 readr_1.3.1                 tidyr_1.0.2                 tibble_2.1.3                tidyverse_1.3.0             cowplot_1.0.0               ggmap_3.0.0                 viridis_0.5.1               viridisLite_0.3.0           kBET_0.99.6                 purrr_0.3.3                 ggpubr_0.2.5                magrittr_1.5                Seurat_3.1.4                scater_1.14.6               ggplot2_3.3.0               SingleCellExperiment_1.8.0  SummarizedExperiment_1.16.1 DelayedArray_0.12.2         BiocParallel_1.20.1         matrixStats_0.55.0          Biobase_2.46.0              GenomicRanges_1.38.0        GenomeInfoDb_1.22.0         IRanges_2.20.2              S4Vectors_0.24.3            BiocGenerics_0.32.0         BiocStyle_2.14.4           
## 
## loaded via a namespace (and not attached):
##   [1] reticulate_1.14          tidyselect_1.0.0         htmlwidgets_1.5.1        grid_3.6.1               Rtsne_0.15               munsell_0.5.0            codetools_0.2-16         mutoss_0.1-12            ica_1.0-2                DT_0.12                  future_1.16.0            withr_2.1.2              colorspace_1.4-1         knitr_1.28               rstudioapi_0.11          ROCR_1.0-7               ggsignif_0.6.0           gbRd_0.4-11              listenv_0.8.0            Rdpack_0.11-1            labeling_0.3             RgoogleMaps_1.4.5.3      GenomeInfoDbData_1.2.2   mnormt_1.5-6             farver_2.0.3             vctrs_0.2.3              generics_0.0.2           TH.data_1.0-10           xfun_0.12                R6_2.4.1                 ggbeeswarm_0.6.0         rsvd_1.0.3               bitops_1.0-6             assertthat_0.2.1         promises_1.1.0           scales_1.1.0             multcomp_1.4-12          beeswarm_0.2.3           gtable_0.3.0             npsurv_0.4-0             globals_0.12.5           sandwich_2.5-1           rlang_0.4.5              splines_3.6.1            lazyeval_0.2.2           broom_0.5.5              BiocManager_1.30.10     
##  [48] yaml_2.2.1               reshape2_1.4.3           modelr_0.1.6             crosstalk_1.0.0          backports_1.1.5          httpuv_1.5.2             tools_3.6.1              bookdown_0.18            gplots_3.0.3             RColorBrewer_1.1-2       ggridges_0.5.2           TFisher_0.2.0            Rcpp_1.0.3               plyr_1.8.6               zlibbioc_1.32.0          RCurl_1.98-1.1           pbapply_1.4-2            zoo_1.8-7                haven_2.2.0              ggrepel_0.8.1            cluster_2.1.0            fs_1.3.2                 data.table_1.12.8        RSpectra_0.16-0          magick_2.3               lmtest_0.9-37            reprex_0.3.0             RANN_2.6.1               mvtnorm_1.1-0            fitdistrplus_1.0-14      xtable_1.8-4             mime_0.9                 hms_0.5.3                patchwork_1.0.0          lsei_1.2-0               evaluate_0.14            jpeg_0.1-8.1             readxl_1.3.1             gridExtra_2.3            compiler_3.6.1           KernSmooth_2.23-16       crayon_1.3.4             htmltools_0.4.0          later_1.0.0              RcppParallel_4.4.4       lubridate_1.7.4          DBI_1.1.0               
##  [95] dbplyr_1.4.2             MASS_7.3-51.5            Matrix_1.2-18            cli_2.0.2                gdata_2.18.0             metap_1.3                igraph_1.2.4.2           pkgconfig_2.0.3          sn_1.5-5                 numDeriv_2016.8-1.1      sp_1.4-1                 plotly_4.9.2             xml2_1.2.2               vipor_0.4.5              multtest_2.42.0          XVector_0.26.0           bibtex_0.4.2.2           rvest_0.3.5              digest_0.6.25            sctransform_0.2.1        RcppAnnoy_0.0.15         tsne_0.1-3               rmarkdown_2.1            cellranger_1.1.0         leiden_0.3.3             uwot_0.1.5               DelayedMatrixStats_1.8.0 shiny_1.4.0              gtools_3.8.1             rjson_0.2.20             lifecycle_0.1.0          nlme_3.1-145             jsonlite_1.6.1           BiocNeighbors_1.4.2      fansi_0.4.1              pillar_1.4.3             lattice_0.20-40          fastmap_1.0.1            httr_1.4.1               plotrix_3.7-7            survival_3.1-8           glue_1.3.1               FNN_1.1.3                png_0.1-7                stringi_1.4.6            BiocSingular_1.2.2       caTools_1.18.0          
## [142] irlba_2.3.3              future.apply_1.4.0       ape_5.3